﻿using System;
using System.Diagnostics;
using Nature.Numerics;
using Nature.Numerics.Basic;

namespace Nature
{
    /// <summary>
    /// Polytropic Gas Ratios
    /// </summary>
    [DebuggerDisplay("compression: {DensityRatio}, T-ratio: {TemperatureRatio}, P-ratio: {PressureRatio}")]
    public struct PolytropicGasRatios : IEquatable<PolytropicGasRatios>
    {
        /// <summary>
        /// P / P0
        /// </summary>
        public double PressureRatio { get; set; }

        /// <summary>
        /// T / T0
        /// </summary>
        public double TemperatureRatio { get; set; }

        /// <summary>
        /// V0 / V (same as the compression ratio)
        /// </summary>
        public double DensityRatio { get; set; }


        public bool Equlas(PolytropicGasRatios other, DoubleRelativeTolerance reltol)
        {
            return
                reltol.AreSame(this.PressureRatio, other.PressureRatio) &&
                reltol.AreSame(this.TemperatureRatio, other.TemperatureRatio) &&
                reltol.AreSame(this.DensityRatio, other.DensityRatio);
        }

        #region Methods Calculating the Mach Number

        public static MachNumber CalcMachNumber(CompressionRatio cr, double gamma)
        {
            double m2 = 2.0 * cr
                / ((gamma + 1.0) - cr * (gamma - 1.0));
            return new MachNumber(Math.Sqrt(m2));
        }

        public static MachNumber CalcMachNumber(Temperature t0, Temperature t1, double gamma)
        {
            double delta = Math.Max(t0.ToDouble(), t1.ToDouble()) / Math.Min(t0.ToDouble(), t1.ToDouble());

            Newton solver = Newton.Create(m =>
            {
                double tmp = 2.0 * gamma * m * m - (gamma - 1.0);
                tmp *= (gamma - 1.0) * m * m + 2.0;
                tmp /= (gamma + 1.0) * (gamma + 1.0) * m * m;
                return delta - tmp;
            });
            solver.RelTol = 1.0e-12;
            solver.AbsTol = 1.0e-12;
            double result = solver.Solve(2.0);
            return new MachNumber(result);
        }

        public static MachNumber CalcMachNumber(Pressure p0, Pressure p1, double gamma)
        {
            double delta = Math.Max(p0.ToDouble(), p1.ToDouble()) / Math.Min(p0.ToDouble(), p1.ToDouble());
            double result = (delta + (gamma - 1.0) / (gamma + 1.0))
                * (gamma + 1.0) / 2.0 / gamma;
            result = Math.Sqrt(result);
            return new MachNumber(result);
        } 
        #endregion

        #region CaclShockWaveRatios Methods

        public static PolytropicGasRatios CaclShockWaveRatios(MachNumber m, double gamma)
        {
            PolytropicGasRatios ratios = new PolytropicGasRatios();
            double m2 = m * m;

            ratios.PressureRatio = (2.0 * gamma / (gamma + 1.0) * m2) - (gamma - 1.0) / (gamma + 1.0);

            ratios.TemperatureRatio =
                (2.0 * gamma * m2 - (gamma - 1.0))
                * ((gamma - 1) * m2 + 2.0)
                / ((gamma + 1.0) * (gamma + 1.0) * m2);

            ratios.DensityRatio = ((gamma + 1.0) * m * m) / ((gamma - 1.0) * m * m + 2.0);

            return ratios;
        }

        public static PolytropicGasRatios CaclShockWaveRatios(CompressionRatio cr, double gamma)
        {
            return CaclShockWaveRatios(CalcMachNumber(cr, gamma), gamma);
        }

        public static PolytropicGasRatios CaclShockWaveRatios(Temperature t0, Temperature t1, double gamma)
        {
            return CaclShockWaveRatios(CalcMachNumber(t0, t1, gamma), gamma);
        }

        public static PolytropicGasRatios CaclShockWaveRatios(Pressure p0, Pressure p1, double gamma)
        {
            return CaclShockWaveRatios(CalcMachNumber(p0, p1, gamma), gamma);
        } 
        #endregion

        #region CalcIsentropicRatios Methods

        public static PolytropicGasRatios CalcIsentropicRatios(CompressionRatio cr, double gamma)
        {
            PolytropicGasRatios ratios = new PolytropicGasRatios();
            ratios.PressureRatio = Math.Pow(cr, gamma);
            ratios.TemperatureRatio = Math.Pow(cr, gamma - 1.0);
            ratios.DensityRatio = cr;
            return ratios;
        }

        public static PolytropicGasRatios CalcIsentropicRatios(Temperature t0, Temperature t1, double gamma)
        {
            PolytropicGasRatios ratios = new PolytropicGasRatios();
            ratios.TemperatureRatio = t1.ToDouble() / t0.ToDouble();
            ratios.DensityRatio = Math.Pow(ratios.TemperatureRatio, 1.0 / (gamma - 1.0));
            ratios.PressureRatio = ratios.TemperatureRatio * ratios.DensityRatio;
            return ratios;
        }


        public static PolytropicGasRatios CalcIsentropicRatios(Pressure p0, Pressure p1, double gamma)
        {
            PolytropicGasRatios ratios = new PolytropicGasRatios();
            ratios.PressureRatio = p1.ToDouble() / p0.ToDouble();
            ratios.DensityRatio = Math.Pow(ratios.PressureRatio, 1.0 / gamma);
            ratios.TemperatureRatio = ratios.PressureRatio / ratios.DensityRatio;
            return ratios;
        } 
        #endregion

        #region IEquatable<PolytropicGasRatios> Members

        public bool Equals(PolytropicGasRatios other)
        {
            return
                this.PressureRatio == other.PressureRatio &&
                this.DensityRatio == other.DensityRatio &&
                this.TemperatureRatio == other.TemperatureRatio;
        }

        #endregion
    }
}
